% Runs simluation of Big G model
% This version: December 2023
% Generates figures 5 and 6
% Calls mod file: modelG.mod

clear all; close all; clc

addpath c:\dynare\5.4\matlab

variables = char('federalpurchases','pi_c','i','y','g_output');
var_names = {'Federal purchases (percent)','Inflation (percent)','Interest Rate (p.p.)','Output (percent)','Spendingshare'};



n_var = size(variables,1);


timemax = 21;
time    = 0:timemax-1;

dynare modelG noclearall           


options_.nomoments=1;
options_.nofunctions=1;
options_.irf=timemax;
options_.nograph=1;
options_.noprint=1;
IRF1=NaN(n_var,timemax);
IRF2=NaN(n_var,timemax);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% baseline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set_param_value('gy',.187);    % output share of total gov spending
set_param_value('rho_g1',.65); % persistence shock 1
set_param_value('rho_g2',.73); % persistence shock 2
set_param_value('std1',13.5);    % std shock 1 (percent)
set_param_value('std2',13.1);    % std shock 2 (percent)
set_param_value('alpha1',.78); % calvo sec 1
set_param_value('alpha2',.89); % calvo sec 2
set_param_value('gamma',1-.76);  % fraction of total spending going to sector 1 (interpolated form federal purchases)
set_param_value('n',0.65);     % size sector 1 (value added)


[info, oo_] = stoch_simul(M_, options_, oo_, var_list_);
if info==0
    for varindex = 1:n_var
        IRF1_base(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g1']);
        IRF2_base(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g2']);
    end
else
    fprintf('Could not solve model for grid point %5.4f',ngrid(ijk))
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% symmetric model 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


averagecalvo = .5*.78+.5*.89;
set_param_value('gamma',.5);
set_param_value('n',.5);
set_param_value('alpha1',averagecalvo);
set_param_value('alpha2',averagecalvo);


[info, oo_] = stoch_simul(M_, options_, oo_, var_list_);
if info==0
    for varindex = 1:n_var
        IRF1_sym(varindex,:) = oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g1']);
        IRF2_sym(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g2']);
    end
else
    fprintf('Could not solve model for grid point %5.4f',ngrid(ijk))
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% figure 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
for ijk = 1:4
    subplot(4,2,(ijk-1)*2 +1)
    plot(time,IRF1_base(ijk,:),'b',time,IRF1_sym(ijk,:),'r--','LineWidth',1.5); hold on;
    plot(time,zeros(1,timemax),'k--')
    title(deblank(var_names(ijk)) )
    ylabel(deblank(var_names(ijk)) )
  %   ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    if ijk == 1
        legend('Baseline','Symmetric model')
        title('Shock sector 1')
    end


    subplot(4,2,(ijk-1)*2+2)
    plot(time,IRF2_base(ijk,:),'b',time,IRF2_sym(ijk,:),'r--','LineWidth',1.5); hold on;
    plot(time,zeros(1,timemax),'k--')
   %  ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    if ijk == 1
        title('Shock sector 2')
    end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% figure 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for ijk  =1:timemax
  discountvector(ijk) = 0.997^(ijk-1); 
end
discountvector = discountvector';

for ijk = 1:timemax
    pvmult1(ijk) = IRF1_base(find(strcmp(var_names,'Output (percent)')),1:ijk)*discountvector(1:ijk)./(IRF1_base(find(strcmp(var_names,'Spendingshare')),1:ijk)*discountvector(1:ijk));
    pvmult2(ijk) = IRF2_base(find(strcmp(var_names,'Output (percent)')),1:ijk)*discountvector(1:ijk)./(IRF2_base(find(strcmp(var_names,'Spendingshare')),1:ijk)*discountvector(1:ijk));
end


figure
subplot(1,2,1)
plot(time,pvmult1,'b-',time,pvmult2,'r:','LineWidth',2);
legend('Sec 1 shock','Sec 2 shock','Interpreter','latex'); legend('BoxOff')
%ylim([0.35 .8])
%xlim([0 24])
grid on
xlabel('Months','Interpreter','latex')
set(gca,'TickLabelInterpreter','latex')
set(gca,'FontSize',30)
fig = gcf;
set(fig,'PaperOrientation','landscape');
%print(gcf, ['../../figures/pvmult.pdf'], '-bestfit', '-dpdf') 



subplot(1,2,2)
plot(time,IRF1_base(find(strcmp(var_names,'Spendingshare')),:),'k',time,IRF2_base(find(strcmp(var_names,'Spendingshare')),:),'k:',time,IRF1_base(find(strcmp(var_names,'Output (percent)')),:),'b-',time,IRF2_base(find(strcmp(var_names,'Output (percent)')),:),'r:','LineWidth',2); hold on
legend('Gov.spending Sec 1 shock','Gov.spending Sec 2 shock','Output Sec 1 shock','Output Sec 2 shock','Interpreter','latex'); legend('BoxOff')
ylabel('Percent of output','Interpreter','latex','FontWeight','normal');%,'FontSize',12,'FontWeight','normal');
%ylim([0 .35])
%xlim([0 24])
 grid on
xlabel('Months','Interpreter','latex')
set(gca,'TickLabelInterpreter','latex')
set(gca,'FontSize',30)
fig = gcf;
set(fig,'PaperOrientation','landscape');
%print(gcf, ['../../figures/pvmult_components.pdf'], '-bestfit', '-dpdf') 


%close all